function [x, y, z] = lla2ned(lat, lon, alt, lat0, lon0, alt0)
% LLA to NED conversion using WGS84 ellipsoid
% Inputs:
%   lat, lon, alt: latitude, longitude, altitude (in degrees and meters)
%   lat0, lon0, alt0: reference point coordinates (in degrees and meters)
% Outputs:
%   x, y, z: North, East, Down coordinates (in meters)
% 
% Example usage:
% [x, y, z] = lla2ned(37.7749, -122.4194, 20, 37.7748, -122.4193, 10);

% WGS84 ellipsoid parameters
a = 6378137.0; % semimajor axis (m)
b = 6356752.314245; % semiminor axis (m)
f = (a - b) / a; % flattening

% Convert reference point to ECEF
N0 = a / sqrt(1 - f^2 * sind(lat0)^2);
X0 = (N0 + alt0) * cosd(lat0) * cosd(lon0);
Y0 = (N0 + alt0) * cosd(lat0) * sind(lon0);
Z0 = (N0 * (1 - f^2) + alt0) * sind(lat0);

% Convert LLA to ECEF
N = a / sqrt(1 - f^2 * sind(lat)^2);
X = (N + alt) * cosd(lat) * cosd(lon);
Y = (N + alt) * cosd(lat) * sind(lon);
Z = (N * (1 - f^2) + alt) * sind(lat);

% Convert ECEF to NED
cosLat0 = cosd(lat0);
sinLat0 = sind(lat0);
cosLon0 = cosd(lon0);
sinLon0 = sind(lon0);
R_ECEFtoNE = [-sinLat0*cosLon0, -sinLat0*sinLon0, cosLat0; ...
              -sinLon0,           cosLon0,           0; ...
              -cosLat0*cosLon0, -cosLat0*sinLon0, -sinLat0];
D_EtoNED = [0, 1, 0; ...
            1, 0, 0; ...
            0, 0,-1];
delta_ECEF = [X - X0; Y - Y0; Z - Z0];
delta = D_EtoNED * R_ECEFtoNE * delta_ECEF;
x = delta(1);
y = delta(2);
z = delta(3);
end